{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Velocity Data Download Script for Validation\n",
    "#### Working Script to pull large velocity datasets and later compare to ICESAT2-derived velocities\n",
    "\n",
    "ICESat-2 hackweek  \n",
    "June 15, 2020  \n",
    "Lynn Kaluzienski"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os,re,h5py\n",
    "import requests\n",
    "import zipfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import pyproj\n",
    "import scipy, sys, os, pyproj, glob\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import Point, Polygon\n",
    "import pointCollection as pc\n",
    "import math\n",
    "\n",
    "# Import some of the scripts that we have written\n",
    "import sys\n",
    "sys.path.append(\"/home/jovyan/surface_velocity/scripts\")\n",
    "from loading_scripts import atl06_to_dict\n",
    "\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import MEaSUREs Velocity Datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-86 -81 -81 -86 -86]\n",
      "[-55 -55 -65 -65 -55]\n"
     ]
    }
   ],
   "source": [
    "#From Ben Smith's code loading in .tif file, running into issues likely with directories\n",
    "data_root='/srv/shared/surface_velocity/FIS_Velocity/'\n",
    "#spatial_extent = np.array([-102, -76, -98, -74.5])\n",
    "spatial_extent = np.array([-65, -86, -55, -81])\n",
    "lat=spatial_extent[[1, 3, 3, 1, 1]]\n",
    "lon=spatial_extent[[2, 2, 0, 0, 2]]\n",
    "print(lat)\n",
    "print(lon)\n",
    "# project the coordinates to Antarctic polar stereographic\n",
    "xy=np.array(pyproj.Proj(3031)(lon, lat))\n",
    "# get the bounds of the projected coordinates \n",
    "XR=[np.nanmin(xy[0,:]), np.nanmax(xy[0,:])]\n",
    "YR=[np.nanmin(xy[1,:]), np.nanmax(xy[1,:])]\n",
    "#Originally tried to load data from a local directory, should change to shared directory\n",
    "Measures_vx=pc.grid.data().from_geotif(os.path.join(data_root,'Measures2_FIS_Vx.tif'), bounds=[XR, YR])\n",
    "Measures_vy=pc.grid.data().from_geotif(os.path.join(data_root,'Measures2_FIS_Vy.tif'), bounds=[XR, YR])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load a line and plot\n",
    "\n",
    "data_root='/srv/shared/surface_velocity/'\n",
    "field_dict={None:['delta_time','latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "\n",
    "rgt = \"0848\"\n",
    "cycle=\"03\"\n",
    "filename = glob.glob(os.path.join(data_root, 'FIS_ATL06', f'*ATL06_*_{rgt}{cycle}*_003*.h5'))[0]\n",
    "\n",
    "D=atl06_to_dict(filename,'/gt2l', field_dict=field_dict, index=None, epsg=3031)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Plot Original MEaSUREs Maps of Velocity Components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show the velocity maps:\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(121)\n",
    "Measures_vx.show(cmap='viridis', clim=[-100,100])\n",
    "plt.plot(xy[0,:], xy[1,:],'k')\n",
    "plt.title('Measures X-Velocity')\n",
    "plt.plot(D['x'],D['y'],'r')\n",
    "plt.gca().set_aspect('equal')\n",
    "\n",
    "plt.subplot(122)\n",
    "Measures_vy.show(cmap='viridis', clim=[-100,100])\n",
    "plt.plot(xy[0,:], xy[1,:],'k')\n",
    "plt.title('Measures Y-Velocity')\n",
    "plt.plot(D['x'],D['y'],'r')\n",
    "plt.gca().set_aspect('equal')\n",
    "\n",
    "#plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolate the Measures velocities along track and rotate Velocities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:16: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  app.launch_new_instance()\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8c3157768c374f698621c8f861276c94",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "00f5c8950f9448f88450ee55c09a8eba",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vx = Measures_vx.interp(D['x'],D['y'])\n",
    "vy = Measures_vy.interp(D['x'],D['y'])\n",
    "\n",
    "#Solve for angle to rotate Vy to be along track and Vx to be across track\n",
    "\n",
    "xL=abs((D['x'][0])-(D['x'][1]))\n",
    "yL=abs((D['y'][0])-(D['y'][1]))\n",
    "import math\n",
    "theta_rad=math.atan(xL/yL)\n",
    "#theta_deg=theta_rad*180/math.pi\n",
    "v_along=vy/math.cos(theta_rad)\n",
    "v_across=vx/math.cos(theta_rad)\n",
    "\n",
    "\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(221)\n",
    "plt.plot(D['x_atc'],vx)\n",
    "plt.title('MEaSUREs Vx')\n",
    "plt.subplot(222)\n",
    "plt.plot(D['x_atc'],vy)\n",
    "plt.title('MEaSUREs Vy')\n",
    "plt.subplot(223)\n",
    "plt.plot(D['x_atc'],v_along)\n",
    "plt.title('V_along')\n",
    "plt.subplot(224)\n",
    "plt.plot(D['x_atc'],v_across)\n",
    "plt.title('V_across')\n",
    "plt.tight_layout()\n",
    "\n",
    "#plt.figure(figsize=(8,4))\n",
    "#plt.subplot(121)\n",
    "#plt.plot(D['x_atc'],v_along)\n",
    "#plt.title('V_along')\n",
    "#plt.subplot(122)\n",
    "#plt.plot(D['x_atc'],v_across)\n",
    "#plt.title('V_across')\n",
    "\n",
    "# Double check Velocities, v_along should be similar to vy but should major differences where flow angle changes\n",
    "Vdiff=vy-v_along\n",
    "\n",
    "\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.plot(D['x_atc'],Vdiff)\n",
    "plt.title('Difference between Vy and Valong')\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "\n",
    "#load rotated velocities into dictionary\n",
    "#D['v_along']=v_along"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
